function [dat,H,coefmat,realdat] = genInterval(data,datareal,mes,anio,type)
%H={};

%x=[1.7 1.285 3.708 6.001 7.569 10.170 8.777 13.756 18.681 20.733 26.088 30.880 37.479 46.230]'
x=load(datareal);
realdat=x;
%dat=[1.7 1.285 NaN NaN NaN 10.170 NaN 13.756 18.681 NaN 26.088 NaN NaN 46.230]'
dat=load(data)';
if isnan(dat(1,1))
    dat(1,1)=x(1,1);
    n=randi([1 size(dat,1)],1,1);
  while isnan(dat(n,1))
      n=randi([1 size(dat,1)],1,1);
  end
    dat(n,1)=NaN;
end

if isnan(dat(size(dat,1),1))
    dat(size(dat,1),1)=0;
end

%dat;
%%%
%a=load(data);
%an=a(:,1);
%ind=find(an==anio);
%bn=a(ind(1,1):ind(size(ind,1),1),2);
%indb=find(bn==mes);
%dat=a(ind(1,1)+indb(1,1)-1:ind(1,1)+indb(1,1)+size(indb,1)-2,type);
%as=dat;% data original
%indanio=ind;
%indmes=indb;
%%%
%generando data aleatoria
%%%
%rd=randi([1 size(dat,1)],9,1)
%for i=1:size(rd,1)
%    dat(rd(i,1))=NaN;
%end
%%%
%%generando intervalos
%%%
num=NaN;
count = 1;
count2=1;

while count<=size(dat,1)
    if ~isnan(dat(count,1)) 
        if isnan(num)
        num=count;
       else 
           H{count2}=[num,count,count-num];
            num=count;
            count2= count2+1;
        end 
    end
    count=count+1;
    
    
end

%So y Sn se consideran
coefmat=zeros(size(H,2)-1,size(H,2)-1);
coefmat(1,1)=2*(H{1}(3)+H{2}(3));
coefmat(1,2)=H{2}(3);


coefmat(size(coefmat,1),size(coefmat,2))=2*(H{size(H,2)-1}(3)+H{size(H,2)-2}(3));
coefmat(size(coefmat,1),size(coefmat,2)-1)=H{size(H,2)-1}(3);

countcoef=1;
for r=2:size(coefmat,1)-1   
        coefmat(r,countcoef)=H{r}(3);
        coefmat(r,countcoef+1)=2*(H{r}(3)+H{r+1}(3));
        coefmat(r,countcoef+2)=H{r+1}(3);
        countcoef=countcoef+1;
end

for l=2:size(H,2)          
    ymat(l-1,1)=(dat(H{l}(2),1)-dat(H{l}(1),1))/(H{l}(3)) - (dat(H{l}(1),1)-dat(H{l-1}(1),1))/(H{l-1}(3)) ;
end

    ymat=6*ymat;
   %%sn es la matris de coeicientes S1 . .. Sn  
    sn=inv(coefmat)*ymat; %S1... Sn
    
    sn([size(sn,1)+1],1) = 0;
    %sn=sn';
    for jj=1:size(coefmat,1)
        S(jj,1)=(sn(jj+1,1)-sn(jj,1))/(6*H{jj+1}(3));
        S(jj,2)=sn(jj,1)/2;
        S(jj,3)=(dat(H{jj+1}(2))-dat(H{jj}(2)))/H{jj+1}(3) - (2*H{jj+1}(3)*sn(jj,1)+H{jj+1}(3)*sn(jj+1,1))/6;
        S(jj,4)=dat(H{jj}(2));
    end
    dat;
    datcom=dat;
    
    
    
    
    
    countt=1;
    for ii=1:size(H,2)
        
        if H{ii}(3)>1 
            for ji=H{ii}(1):H{ii}(2)-2
            %fprintf('entrea a %d',ji);
             if ii ~=1 
             datcom(ji+1,1)=(S(ii-1,1)*(((ji+1)-H{ii}(1))^3))+(S(ii-1,2)*(((ji+1)-(H{ii}(1)))^2))+(S(ii-1,3)*((ji+1)-H{ii}(1)))+S(ii-1,4);
             else
                datcom(ji+1,1)=(S(ii,1)*(((ji+1)-H{ii}(1))^3))+(S(ii,2)*(((ji+1)-(H{ii}(1)))^2))+(S(ii,3)*((ji+1)-H{ii}(1)))+S(ii,4);
             end
               countt=countt+1;
            end
        end        
    end
        
        %RMSE
        error=sqrt(sum((datcom-x).^2)/(size(dat,1)));
        
        error
        %MSE
        ms=mse(datcom-x)
        xx=[1:size(dat,1)];
      h1=plot(xx,x,'-rs');hold on;
       h2=plot(xx,datcom,'-bs');
     legend([h1 h2],{'real','cubicspline'});
       
        
        
        
        

end